Matrix Operations

Please refer Computational Physics-Numerical-Solution of Linear Equations
LUP Decomposition
PA=LU about Ax=b →PAx=Pb →LUx=Pb →Ly=Pb with y=Ux (전진대입) →Ux=y (후진 대입) 전진대입 Θ(n2) yi=bπ[i]−Σi−1j=1lijyj 후진대입 Θ(n2) xi={yi−Σnj=i+1uijxj}/uii
#include <iostream>
#include <vector>
using namespace std;
typedef vector<float> vec1d;
typedef vector<vector<float>> vec2d;
vec1d LUP_solve(vec2d L, vec2d U, vec1d pi, vec1d b){
int n=L.size();
vec1d x(n);
vec1d y(n);
for(int i=0; i<n; ++i){
float ly=0;
for(int j=0; j<i; ++j){
ly+=L[i][j]*y[j];
}
y[i]=b[pi[i]]-ly;
}
for(int i=n-1; i>=0; --i){
float ux=0;
for(int j=0; j<n; ++j){
ux+=U[i][j]*x[j];
}
x[i]=(y[i]-ux)/U[i][i];
}
return x;
}
void LU_Decomposition(vec2d A, vec2d& L, vec2d& U){
int n=A.size();
for(int i=0; i<n; ++i){
for(int j=0; j<n; ++j){
if(i>j) U[i][j]=0;
else if(i<j) L[i][j]=0;
else L[i][i]=1;
}
}
for(int k=0; k<n; ++k){
U[k][k]=A[k][k];
for(int i=k+1; i<n; ++i){
L[i][k]=A[i][k]/A[k][k];
U[k][i]=A[k][i];
}
for(int i=k+1; i<n; ++i){
for(int j=k+1; j<n; ++j){
A[i][j]=A[i][j]-L[i][k]-U[k][j];
}
}
}
}
int main(void){
vec2d A={
{1, 2, 0},
{3, 4, 4},
{5, 6, 3}
};
vec1d b={3, 7, 8};
vec2d L(3, vec1d(3));
vec2d U(3, vec1d(3));
LU_Decomposition(A, L, U);
vec1d P={2, 0, 1};
vec1d x=LUP_solve(L, U, P, b);
for(int i=0; i<A.size(); ++i){
for(int j=0; j<A.size(); ++j) std::cout<<L[i][j]<<' ';
std::cout<<'\n';
}
for(int i=0; i<A.size(); ++i){
for(int j=0; j<A.size(); ++j) std::cout<<U[i][j]<<' ';
std::cout<<'\n';
}
for(float xx: x) std::cout<<xx<<' ';
std::cout<<std::endl;
return 0;
}
LU 분해는 가우스 소거법으로 얻을 수 있음
#include <iostream>
#include <vector>
using namespace std;
typedef vector<float> vec1d;
typedef vector<vector<float>> vec2d;
void LU_Decomposition(vec2d A, vec2d& L, vec2d& U){
int n=A.size();
for(int i=0; i<n; ++i){
for(int j=0; j<n; ++j){
if(i>j) U[i][j]=0;
else if(i<j) L[i][j]=0;
else L[i][i]=1;
}
}
for(int k=0; k<n; ++k){
U[k][k]=A[k][k];
for(int i=k+1; i<n; ++i){
L[i][k]=A[i][k]/A[k][k];
U[k][i]=A[k][i];
}
for(int i=k+1; i<n; ++i){
for(int j=k+1; j<n; ++j){
A[i][j]=A[i][j]-L[i][k]*U[k][j];
}
}
}
}
int main(void){
vec2d A={
{2, 3, 1, 5},
{6, 13, 5, 19},
{2, 19, 10, 23},
{4, 10, 11, 31}
};
vec2d L(4, vec1d(4));
vec2d U(4, vec1d(4));
LU_Decomposition(A, L, U);
for(int i=0; i<A.size(); ++i){
for(int j=0; j<A.size(); ++j) std::cout<<L[i][j]<<' ';
std::cout<<'\n';
}
std::cout<<std::endl;
for(int i=0; i<A.size(); ++i){
for(int j=0; j<A.size(); ++j) std::cout<<U[i][j]<<' ';
std::cout<<'\n';
}
return 0;
}
대각 원소가 0이 될 때, 0으로 나누어줘야 하는 경우가 발생할 수 있다.
이 때, 피벗 P를 이용하면 항상 분해 가능하다.
#include <iostream>
#include <vector>
#include <algorithm>
#define swap(i, j){ int tmp=i; i=j; j=tmp; }
using namespace std;
typedef vector<float> vec1d;
typedef vector<vector<float>> vec2d;
void LUP_Decomposition(vec2d A, vec1d& P){
int n=A.size();
for(int i=0; i<n; ++i){
P[i]=i;
}
for(int k=0; k<n; ++k){
int kp=0;
float p=0;
for(int i=k; i<n; ++i){
if(std::abs(A[i][k])>p){
p=std::abs(A[i][k]);
kp=i;
}
}
if(p==0){
perror("Singular Matrix");
exit(1);
}
swap(P[k], P[kp]);
for(int i=0; i<n; ++i){
swap(A[k][i], A[kp][i]);
}
for(int i=k+1; i<n; ++i){
A[i][k]=A[i][k]/A[k][k];
for(int j=k+1; j<n; ++j){
A[i][j]=A[i][j]-A[i][k]*A[k][j];
}
}
}
}
int main(void){
vec2d A={
{2, 0, 2, .6},
{3, 3, 4, -2},
{5, 5, 4, 2},
{-1, -2, 3.4, -1}
};
vec1d P(4);
LUP_Decomposition(A, P);
for(float p: P) std::cout<<p<<' ';
std::cout<<std::endl;
return 0;
}